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Abstract 

We study the long-time behavior of the scaled walker (particle) position associated with de- 
coupled continuous-time random walk which is characterized by superheavy-tailed distribution of 
waiting times and asymmetric heavy-tailed distribution of jump lengths. Both the scaling function 
and the corresponding limiting probability density are determined for all admissible values of tail 
indexes describing the jump distribution. To analytically investigate the limiting density function, 
we derive a number of different representations of this function and, by this way, establish its main 
properties. We also develop an efficient numerical method for computing the limiting probability 
density and compare our analytical and numerical results. 
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I. INTRODUCTION 



The continuous-time random walks (CTRWs), i.e., cumulative jump processes which are 
characterized by a joint probability density of waiting time and jump length, play a signifi- 
cant role in many areas of science. The reason is that a large variety of physical, biological 
and other systems are often modeled by two random variables that can be interpreted as the 
waiting time between successive transitions (jumps) of the system into new states and the 
transition measure (jump length). For example, the CTRW model can be used to describe 
anomalous diffusion and transport in disordered media [1-4], human mobility [5, 6], financial 
[7-10] and seismic [11, 12] data. 

According to the theory of CTRWs [13] (see also Refs. [2, 3]), the probability density 
P(x,t) of the particle position X(t) depends only on the joint probability density of wait- 
ing time and jump length. Unfortunately, even in the simplest (decoupled) case when the 
joint density is a product of waiting-time density p(r) and jump-length density w(^), the 
representation of P(x,t) in terms of special functions is known in a few cases [14-16]. In 
contrast, the long-time behavior of P(x, t) that determines the diffusion and transport prop- 
erties of walking particles is studied analytically in much more detail [17-21]. Specifically, 
the asymptotic behavior of the scaled particle position Y(t) = a{t)X{t) is investigated in 
Refs. [20, 21]. In these works, the scaling function a(t) and the distribution function of Y(t) 
at t — > oo are obtained for all waiting-time and jump densities characterized by finite second 
moments or heavy tails. 

A special case of CTRWs with superheavy-tailed distributions of waiting time was first 
considered in [22]. Because all fractional moments of these distributions are infinite, they can 
be used to model extremely anomalous behavior of systems and processes as, e.g., iterated 
maps [23], ultraslow kinetics [24], superslow diffusion [25], and Langevin dynamics [26, 27]. 
For the CTRWs characterized by arbitrary superheavy-tailed distributions of waiting time, 
the scaling function a{t) and the corresponding limiting probability density V(y) of Y(oo) 
have already been found for the jump densities having finite second moments [28] and 
symmetric heavy tails [29]. In the present work, we report a comprehensive theoretical 
and numerical studies of the long-time behavior of the reference CTRWs in the general case 
of asymmetric jump densities characterized by heavy tails. 

The paper is organized as follows. In Sec. II, we describe the CTRW model, formulate 
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the problem of the limiting probability density and list some previously obtained results. In 
Sec. Ill, we find the scaling function a(t) and the limiting probability density V(y) in terms of 
the inverse Fourier transform for all possible cases. A number of alternative representations 
of the limiting density and its main properties are obtained in Sec. IV. Here, we derive 
V{y) in terms of (i) the inverse Mellin transform, (ii) the Laplace transform, (iii) the Fox H 
function, and (iiii) the series expansion. Using the series and Laplace representations of the 
limiting probability density, in Sec. V we determine its short- and long-distance behavior. 
In Sec. VI, we develop a method for the numerical evaluation of V(y) and compare the 
analytical and numerical results. Our findings are summarized in Sec. VII. Finally, a short 
derivation of the fractional equation for V(y) is presented in the Appendix. 



II. MODEL AND BACKGROUND 



One of the main statistical characteristics of a CTRW is the probability density P(x, t) of 
the particle position X(t). This quantity is completely determined by the joint probability 
density t) of waiting times r n (r„ > 0), i.e., times between successive jumps, and 
jump lengths (—00 < < 00) that are assumed to be independent and identically 
distributed random variables with the probability densities p{r) and w(£), respectively. In 
the case of decoupled CTRWs, when the sets {r n } and {£ n } are independent of each other, 
r) = w(£)p(r) and the probability density P(x,t) in Fourier-Laplace space is given by 
the Montroll- Weiss equation [13] 

P ks = l ~ Ps (2.1) 
s(l-p s w k ) 



Here, w k = J c {w(x)} = dxe tkx w(x) with —00 < k < 00 is the Fourier transform of 
Ps — £{p(t)} — /o°° dte~ st p(t) with Res > is the Laplace transform of p(r), and 
P ks = F{£{P(x,t)}}. 

Representing Eq. (2.1) in the form 

P ks = 1 ~ Ps + ^ ~ Ps ^ PsWk (2 2) 

and taking the inverse Fourier transform [defined as J^^ifk} = f(x) = 
(2n)' 1 JZ Q dke-^f k ] of Eq. (2.2), one gets 

P s (x) = 1 -^S(x) + 8-P^Ljr-i!^\ (2.3) 
s s [1- p s w k J 



with S(x) being the Dirac 5 function. Then, applying the inverse Laplace transform [defined 
as C~ 1 {g s ) = g(t) = (27ri) _1 J^-^ dse st g s , c is a real number that exceeds the real parts 
of all singularities of g s ] to Eq. (2.3), for the probability density of the particle position we 
obtain 

P(x,t) = V(t)5(x) + W ^-^ W— g*— U, (2.4) 

I S { 1 - p s W k J J 

where 

V{t) = = J^drp(r) (2.5) 

is the survival probability, i.e., the probability that a walking particle remains at the initial 
state X(0) = up to time t. According to the definition (2.5), this probability satisfies the 
conditions V(t) — > as t — > oo and V(t) — >■ 1 as t — > 0. 

Since there are no boundary conditions keeping a walking particle inside a finite region, 
the condition P(x, t) — > as t — > oo must hold for all x. The fact that P(x, t) vanishes in 
the long-time limit suggests to introduce the scaled particle position Y(t) = a(t)X(t), where 
a(t) is a positive scaling function, and study, instead of P(x,t), the asymptotic behavior of 
the probability density V(y,t) = P(y/a{t),t)/a{t) of Y(t). It is clear that if a(t) at t — > oo 
tends to zero fast enough, then the limiting probability density 

V(y) = lim P ( -4t, t] (2.6) 
yy> t^ooa{t) \a(t) ) v 1 

reduces to the degenerate density V(y) = S(y). In contrast, if a(t) at t — > oo grows or 
tends to zero slowly enough, then V(y,t) vanishes in the long-time limit. Finally, for a 
certain asymptotic behavior of a(t) the limiting probability density becomes nonvanishing 
and nondegenerate. The last property is of particular interest because in this case the pair 
V(y) and a(t) determines the asymptotic behavior of the probability density of the particle 
position, P(x,t) ~ a(t)V(a(t)x) as t — > oo. To avoid any misunderstanding, we note that 
this statement is correct only if V(y) ^ 0; in those regions of y where V(y) = (for details, 
see below) the use of the limiting probability density for determining the long-time behavior 
of P(x, t) becomes impractical. 

The problem of finding the pairs V(y) and a(t) has been solved for all typical distri- 
butions of waiting times and jump lengths characterized by both finite second moments 
and heavy tails [20, 21]. Recently, we have partially solved this problem for a new class of 
CTRWs with superheavy-tailed distributions of waiting times [28, 29]. These distributions 
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are characterized by the following asymptotic behavior of the waiting-time density: 

p(r)~^ (t->oo), (2.7) 
r 

where h(r) is a slowly varying function defined by the condition h(fir) ~ h{r) (r — > oo) 
holding for all /i > 0. Since p(r) is normalized, / °° drpij) = 1, the function /i(r) must 
decrease in such a way that h(r) = o(l/ lnr) as r — >■ oo. The main feature of these densities 
is that their fractional moments J °° drr c p{r) are infinite for all c > 0. It has been shown [27] 



'o 
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-oo 



that if p(r) is superheavy-tailed and u>(£) has a finite second moment / 2 = f^d^w^), 
then 

V(y) = ^^e-MH(l iy ) (2.8) 

and 

(y/2V(t)/l 2 , Zi = 

a(f) ~ <J (2.9) 



as £ — >■ oo. Here, 5 at b and are the Kronecker delta (5 0) & = 1 if a = 6 and 5 a b = 

if a 7^ 6) and the Heaviside unit function [H(x) — 1 if x > and = if a; < 0], 

respectively, and l\ = J™ d££w(£) is the first moment of w(£). We note that if Zi 7^ 0, then 
the limiting probability density is one-sided: V(y) = on that semi-axis of y where l±y < 0. 
If the jump density is symmetric, w(—£) = w(£), and has heavy tails, then 

™(0~^ m^oc), (2.10) 

where u > and a G (0,2] is the tail index. According to [29], in this case the limiting 
probability density can be represented in the forms 

(l-l/a,l/a),(l/2,l/2) 
(0,1), (1 -l/a,l/a), (1/2, 1/2) 



p(v) = -h% 

a 



1 /""cfae-l"!* 



2a 



Wo ' 1 + 2cos(7ra/2)a; a + x 

(iJ 2 ' 3 [-] is a particular case of the Fox function) and the corresponding scaling function is 
given by 

Y(l +a )sin(W2) y (t) y /a < « < 2 

a(t) ~ <( v ™ 7 (2.12) 



itln[l/V(t)] • 



a = 2 



(t — > oo) with r(l + a) being the T function. Interestingly, Eq. (2.11) at a = 2 reduces to 
Eq. (2.8) with l\ = 0. But since in this case l 2 = oo, the scaling function in Eq. (2.12) differs 
from that given in Eq. (2.9). It should also be stressed that if p(r) is superheavy-tailed, 
then the survival probability V(t) is a slowly varying function [27]. Therefore, in accordance 
with Eqs. (2.9) and (2.12), the long-time evolution of P(x,t) occurs very slowly. 

In this paper, we will study analytically the long-time solutions of the CTRWs charac- 
terized by both superheavy-tailed distributions of waiting time, whose asymptotic behavior 
is described by Eqs. (2.7), and heavy-tailed distributions of jump length. The last distribu- 
tions are assumed to be asymmetric and can have one or two heavy tails. Since the limiting 
probability density V(y) under certain conditions is determined by the heaviest tail (see be- 
low), we consider, without loss of generality, the jump densities with two heavy tails, whose 
asymptotic behavior is given by 

(£^±oo) (2.13) 

(u± > 0, a± G (0, 2]). In contrast to Ref. [29], now we are concerned with the effects arising 
from the asymmetry of w(£). In addition, we are going to develop a numerical method for 
the simulation of these CTRWs and apply it to verify the theoretical predictions. 



III. LIMITING PROBABILITY DENSITIES AND CORRESPONDING SCALING 
FUNCTIONS 

A. Inverse Fourier transform representation of V{y) 

Under certain conditions (see below), the long-time behavior of the probability density 
P(x,t) is determined by the small-s behavior of the Laplace transform P s (x). In turn, 
according to Eq. (2.3), the last behavior is governed by the small-s behavior of p s . Taking 
into account that 

POO 

1- Ps = dt(l - e- st )p(t) 
Jo 

dqe' q V(q/s) (3.1) 
and using the fact that the survival probability V(t) is a slowly varying function, we obtain 

l-p,~V(l/s) (3.2) 



f 

Jo 



as s — > 0, and Eq. (2.3) in this limit yields 



P ^Ym Hx) + n^)ii-v(i/ s)] 



s 2ns 

oo „—ikx 



x / dk—— . (3.3) 

V(l/s) + l-w k 



oo 



Applying to P s (x) the Tauberian theorem [30] [it states that if the function v(t) is ulti- 
mately monotonic and v s ~ s _7 L(l/s) (0 < 7 < 00) as s — > 0, then v(t) ~ f 1 ~ 1 L(t)/T{ , y) as 
t — > 00, where L(t) is a slowly varying function at infinity], from Eq. (3.3) one gets 

/oo ^—ikx 
dk-— — (3.4) 
V(t) + l-w k 1 ; 

(t — > 00). With this result, we can represent the limiting probability density (2.6) as the 
inverse Fourier transform 

1 r°° e - iK v 

where 

•M-fe 1 !^- < 3 - 6 > 

We remark that, since f*™ dye~ lKy = 2~k5(k) and $(0) = 0, this probability density is 
properly normalized: dyV(y) = 1. It should be noted also that V(y) satisfies a simple 
space-fractional equation (see the Appendix). 

Below, using Eqs. (3.5) and (3.6), we determine the limiting probability density and 
corresponding scaling function in all possible situations. 



B. Jump densities with l\ 7^ 

If the first moment l\ of the jump density exists and is non-zero, then directly from the 
relation 1 — w k = /^(l — e lkx )w(x) one finds 1 — w k ~ —il\k, and thus Eq. (3.6) reduces to 

$(«) = -i«sgn(Z 1 )^!y^ (3.7) 

[sgn(x) = ±1 if x ^ 0]. Choosing the asymptotic behavior of the scaling function in the 
form 



(t — > oo), Eq. (3.7) yields $(/t) = — msgn(Zi). Then, using this result, from Eq. (3.5) we 
obtain the one-sided exponential density 

V(y) = e-^H(l lV ). (3.9) 

This limiting probability density describes all CTRWs characterized by superheavy-tailed 
distributions of waiting time and jump distributions having non-zero first moments. It should 
be emphasized that a class of these jump distributions contains both the distributions with 
finite second moments, see Eq. (2.8), and the heavy-tailed distributions with a± G (1,2]. 



C. Jump densities with a G (1,2) and Zi = 

If the first moment of u>(£) exists and equals zero, then, to find the the asymptotic 
behavior of 1 — Wk as k — > 0, it is reasonable to use the following exact formula: 

If 00 ( x \ 

1 — Wk = 777 / dx(l — cosa;)w + ( — 
\k\ Jo \\k\J 

+— J dx{x — sin a;)w"^— ^ , (3.10) 

where 

w ± {t) = w{t)±w{-t). (3.11) 

Let us introduce the notation 

a = min{o; + , «_} (3-12) 
and consider the case with a G (1,2). Then, using the standard integrals [31] 



°° 1 - COSX 7T 

dx -T— = —- ■ -r-r-7— 7^ (3.13) 







x i+u 2T(1 + i/) sin(7ri//2) 



(0 < v < 2) and 



x — sin a; 



dX x^ 2T(1 + 1/) cos(ttz//2) (3,14) 



(1 < i/ < 2), it can be shown from Eq. (3.10) that 

1 - w k ~ q\k\ a - isgn(k)r\k\ a (3.15) 

(k — > 0), where 

7T 

9 = orTT^ — T^T 77v:( M +^ + + uS aa _) (3.16) 

2r(l + a) sin(7ra/2) 
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and 

7T 

r = or /, , — \ 7 j7w(u+S aa+ - (3.17) 

21 (1 + a) cos(7ra/2) 

It is worth to emphasize that Eqs. (3.15)-(3.17) hold for all admissible values of the largest 
tail index « max = max{a + ,a_}, i.e., for a max £ [a, 2]. 

Assuming that the long-time behavior of the scaling function is given by 

1/q 

(i — >■ oo), from Eqs. (3.6) and (3.15) one obtains 



' ■' I -7= I (3-18) 



= (cosy? — isgn(«;) sin (3.19) 
and the limiting probability density (3.5) in this case can be represented as follows: 

[V> n J K 1 + 2cosym a + K 2a ' ( ' ' 



Here, 



and 



cos ip = y=L : , sin ip = (3.21) 
A/g 2 + r 2 A/g 2 + r 2 



a/^ 2 + r 2 = 



7T 



2r(l + a) sin(7ra/2)| cos(7ra/2)| 
x (2 cos(ira)u + U-5 aa+ 5 aa ^ 



+ u 2 + 5 aa+ + u 2 _5 aa _) 1/2 - (3.22) 

The limiting probability density (3.20) and the corresponding scaling function (3.18) 
describe both symmetric and asymmetric CTRWs. In particular, if a + = a_ = a G (1,2) 
and u + = U- = u, then <p = 0, 

A /^T7 2 = — ^— — , (3.23) 

vy r(l + a)sin(7ra/2)' v ; 

and Eqs. (3.20) and (3.18) are reduced to Eqs. (2.11) and (2.12) describing the symmetric 
CTRWs. 
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D. Jump densities with a £ (0, 1) 



Since at a £ (0, 1) the first moment of the probability density w(£) does not exist, in this 
case it is convenient to use, instead of Eq. (3.10), the following formula: 



j-^j- J dx(l — cos x)w + ^j^ S j 

r 

kJo 



l-w k - 



dx smxw_yjj^j . (3.24) 

At first sight, there are two different situations when a max £ [ex, 1) and cc max £ [1,2]. How- 
ever, because we need to know only the leading term of the asymptotic expansion of 1 — w k 
as k — > 0, we can restrict ourselves to considering Eq. (3.24) for a max £ [a, 1). In this case, 
using the asymptotic formula (2.13) and the standard integrals (3.13) and 



sin x it 



o dX x^+- - 2T(1 + u) cos(ttz//2) (3 ' 25) 
(0 < v < 1), one can show that Eq. (3.24) at k — > reduces to Eq. (3.15) with the 
parameters q and r given by the same Eqs. (3.16) and (3.17). Since these results hold also 
for a max £ [1,2], it can be concluded that the expressions (3.18) and (3.20) for the scaling 
function a(t) and the limiting probability density V(y) are valid not only for a £ (1, 2) but 
also for a £ (0,1). We note that Eqs. (3.18) and (3.20) at a + = = a £ (0,1) and 
u + = u- = u are reduced to Eqs. (2.12) and (2.11), respectively. 

E. Jump densities with a = 1 

Denoting the first and second terms on the right-hand side of Eq. (3.24) by J\ and J 2 , 
respectively, at a — 1 and k — > we obtain 



7T 

Ji ~ -(u + 5 la+ + u-5 la _)\k\ (3.26) 



and 



J 2 ~i/c / dx^^(u + 5 la+ - u^5i a _) 

Jc\k\ x 



lc\k\ 

~ i(u + 5 la+ - u^5 la _)k\n (3.27) 



where c is a positive constant. If the parameter 

P = u + 5 la+ - u_5 la _ (3.28) 
10 



is not equal zero, the term J\ can be neglected in comparison with J 2 . In this case 

1 

\k 



1 — Wk ~ — ipk In — — (3.29) 



(k ->■ 0) and 



= — ipn lim Q ^ In • 



t^oo a(t) 
Assuming that in the long-time limit the condition 



t^oo V(t) \K\a{t) 
-ipK lim ^\ In (3.30) 



P|a(t) ln4-~1 (3.31) 



V(t) a(t) 
holds, i.e., 

»W - (3 - 32) 

as i — )> oo, one easily finds that <&(/c) = — iKSgn(p) and 

7>(y) = e-MH(py). (3.33) 

The comparison of Eqs. (3.33) and (3.9) shows that the limiting probability density at 
a = 1 and p ^ has the same form as in the case of jump densities with l\ = p. Thus, the 
parameter p plays here the role of the first moment of w(£). We note, however, that this 
analogy is not complete because at a — 1 the first moment l\ does not exist. The difference 
between the scaling functions (3.32) and (3.8), which correspond to the limiting densities 
(3.33) and (3.9), has the same origin. 

In the opposite case, when p = (this takes place only if a + = a_ = 1 and u + = u_ = u), 
the contribution of J\ becomes dominant and thus 

1 — w k ~ iiu\k\ (3.34) 

(k ->■ 0) and 

$(k) = nun lim (3.35) 

t^oo V(t) 

Choosing the long-time behavior of the scaling function in the form 

a(t) ~ (3.36) 

7VU 

from Eq. (3.5) with <&(/c) = |/c| we immediately get 

, N 1 f 00 , cos(yre) 

This result is a particular case of Eq. (2.11). 
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F. Jump densities with a = 2 



Since a± < 2, the condition a = 2 implies that a + = a_ = 2. It is clear that if /i 7^ 0, 
then the limiting probability density is given by Eq. (3.9). In contrast, at l± — from 
Eqs. (3.10) and (3.6) one obtains 

1 — cos x 



POO 

1 — Wk ~ (u + + u_)k 2 I dx 

Jc\k\ 



c\k\ X 3 



l-NJ 



1 1 

-(u+ + w_)£; 2 ln — (3.38) 

2 I A; I 

(/c — )• 0, c is a positive parameter) and 

* W = I(„ ++ „_)^ta^k-L ,3. 39) 

If the asymptotic behavior of the scaling function is governed by the relation 



a (t) ~ 2W — - r. (3.40) 

(t ->■ 00), then = k 2 and Eq. (3.5) yields 

^/ n 1 /"°° , cos(y/t) 1 i„i , , , 

— X = 2 e ~ (3 ' 41 > 

The same result was obtained in Ref. [29] under the condition that heavy-tailed jump 
densities with a = 2 are symmetric. However, since the condition l ± — does not imply 
that w(—£) = w(£), the two-sided exponential density (3.41) corresponds to a more wide 
class jump densities characterized by the conditions a = 2 and l\ = 0. It should be noted 
that the limiting probability density (3.41) as well as the limiting density (3.37) can also 
be obtained from the general representation (3.20) by taking the limits a — > 2 and a — > 1, 
respectively. But since the scaling functions (3.40) and (3.32) do not follow from Eq. (3.18), 
we considered these cases separately. 

Thus, according to the above analysis, the CTRWs with superheavy-tailed distributions 
of waiting time are characterized by two different limiting probability densities. The first 
one, the one-sided exponential density (3.9), corresponds to the jump densities whose first 
moments are non-zero. In addition, the same limiting density appears in a special case of 
heavy-tailed jump densities with p = u + 5\ a+ — uSi a _ 7^ 0. The second one, the limiting 
probability density (3.20), describes the reference CTRWs that are characterized by heavy- 
tailed distributions of jump length. If these distributions are symmetric, then the limiting 
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probability density (3.20) is reduced to the symmetric one (2.11), whose properties is well 
established [29]. As for the non-symmetric case, it has never been studied. However, because 
of the oscillating character of the integrand, the use of the limiting density V(y) in the form 
of Eq. (3.20) [we recall that this form of V(y) corresponds to a e (0, 1) or a 6 (1, 2) 
and l\ = 0] is not always convenient. Therefore, to gain more insight into the analytical 
properties of V(y), next we derive its different representations. 



IV. ALTERNATIVE REPRESENTATIONS OF V(y) 

A. Representation of "P(y) in terms of the inverse Mellin transform 

To derive alternative expressions for the limiting probability density (3.20), we first 
rewrite it in the form 

V(y) = V 1 (y) + sgn(y)V2(y), (4.1) 

where 

= I r^ (l + cosy^)cosM 
Kyj 7T J l + 2cos<pK, a + K 2a V ; 



and 

+ 2 COS ip K a + K 2a 



simp f°° K a sm(\y\K) 

V 2 (y) = — J dK - ; - ; , n ( U) 



are even functions of y. Then we calculate the Mellin transform of the functions V ± (y) = 
Vi(y) ±P 2 (y), which represent V(y) for positive and negative y: V(y)\ y ^o = T' ± {y). Using 
the definition of the Mellin transform of a function f(y), f r = Ai{f(y)} = dyf(y)y r ~ 1 , 
and the relation f r = u r v\- r that holds for the function f(y) = J °° dxu(yx)v(x) [32], we 
obtain 

Vf = M{cosy}F 1 _ r ± M{smy}G^ r . (4.4) 

Here, according to Ref. [33], A^{cosy} = T(r) cos(7rr/2) (0 < Rer < 1), A^{siny} = 
r(r) sin(7rr/2) (-1 < Rer < 1), 

1 f° 1 + cos <py a r 
F 1 _ r = - / dy— — — — —y 

7r J 1 + 2 cos <py a + y 2a 

1 f 00 , 1 + cospy i-r j 

= — / dy -y « 

na J 1 + 2 cos p y + y 2 

cos[</?(l — r)/a] 
a sin[7r(l — r)/a] 



(4.5) 
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(1 — a < Rer < 1), and 



a—r 



tt J Q 1 + 2 cos (p y a + y 2a 

1 — r 

/*OCi — — 

smy f dy y ° 

na J 1 + 2 cos (py + y 2 
sin[</?(l — r)/a] 



asin[7r(l — r)/a] 



(4.6) 



(1 - a < Rer < 1 + a). 

Collecting the above results, from Eq. (4.4) one finds 

v ± = r(r)sin[(7ra/2±y)(l-r)/a] 
a sin[7r(l — r)/a] 

with max(l — a, 0) < Rer < 1. From this, using the definition of the inverse Mellin 
transform, M' 1 ^} = f(y) = (2m)- 1 / c c _J~ drf r y~ r , the relation P(y)|^ = ^(j/) and 
Eq. (4.7), we obtain the limiting probability density (3.20) in terms of the inverse Mellin 
transform 

V (v) = TT- / dr . ) !_ r x a \y\ , (4.8) 



where max(l — a, 0) < c < 1 and 



wot 

<i>(y) = -y + s s n (^V- ( 4 - 9 ) 



It is worth to note that both representations of V(y), (3.20) and (4.8), are valid for all 
values of the lowest tail index a from the interval (0,2]. But since the limiting probability 
densities at a = 1 and a = 2 have already been determined in Sees. Ill E and III F, further 
we examine Eq. (4.8) for a G (0,1) and a G (1,2) only. There are four different cases 
associated with these intervals, which we consider separately below. 



1. a £ (0, 1), a + 7^ ct_ 

In this case, Eq. (3.21) together with Eqs. (3.16), (3.17) and (3.22) yields cosy? = 
cos(7ra/2) and siny? = (5 aa+ — S aa _) sin(7ra/2). From the last two equations it follows 
that (p = (5 aa+ — 5 aa _)Tca/2, and Eq. (4.9) reduces to 

TTOt 

0(y) = [l + sgn(<7y)]— , (4.10) 
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where a = 5 aa+ -5 aa _ = sgn(a_ -a+). Denoting <p sgn (ay) = <t>{y), Eq. (4.10) yields 0+ = na 
and 0_ = 0. The last condition means that V(y) — as ay < 0. In other words, in the case 
when a G (0, 1) and a + ^ the limiting probability density V(y) is one-sided. According 
to Eq. (4.8), it can be represented as 



2m Jc-ioo «Sin(7T — J 



with + = TV a. 



2. a £ (1, 2), a+ / a_ ; h = 

For these conditions, Eq. (3.21) leads to the equations cosy? = — cos(7ra/2) and simp = 
— asm(Tca/2), whose solution is given by (p = — a(ir — na/2). In this case Eq. (4.10) reads 

4>{v) = -y - sgn(ay) [tt - — J , (4.12) 

and so <p + = n(a — 1) and 0_ = 7r. Therefore, using the inverse Mellin transform of the T 
function [33] 

— / drT(r)\y\- r = e- lyl (4.13) 
(c > 0), Eq. (4.8) can be rewritten as 



2 ™ Jc-ioo asm(7r— J 
1 



— r 



+ H(-ay)-e- lvl , (4.14) 



a 



where + = n(a — 1). 

Thus, if a G (1,2), a + ^ a_ and l\ = 0, then, in contrast to the previous case, the 
limiting probability density is two-sided. As is clear from Eq. (4.14), this density exhibits 
an exponential decay at y > if a = a_ or at y < if a = a + . We note also that, to avoid 
the double contribution of the point y — 0, the condition H(ya)\ y= o = H(a) is assumed to 
hold. 
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3. a + = ct_ = a G (0, 1), u + / u_ 



Under these conditions, Eq. (3.21) can be expressed as 

(u + + U-) cos(7ra/2) 



cos (p 



sin tp 



Introducing the notation 
from Eq. (4.15) one obtains 



\/u\ + u 2 _ + 2 cos(7tq;)m + m_ ' 

(u+ — m_) sin(7ra/2) 
a/w+ + u 2 _ + 2 cos(7To;)m + m_ 

M + — M_ 



(4.15) 



(p = sgn(e) arctan 



e tan — 
11 V 2 



)]■ 



(4.16) 



(4.17) 



where arctan(x) denotes the principal value of the inverse tangent function. Finally, repre- 
senting the two- valued function (4.9) as <f>(y) = (f) sgn ( e y), where 



i± = — ± arctan 
± 2 



{ixa>\ 
\e\t^{-) 



(4.18) 



and using Eq. (4.8), we find the following two-sided limiting probability density: 



V{y) = 



H{ey) 
2m 



c+ioo r ( r ) gin 

dr — 



a ) 



a sin 



y\ 



H(-ey) 
2m 



c+ioo 



dr- 



T(r) sin I 



a 



sin^) 



l-r \ 
a J I I — r 

y\ 



(4.19) 



Since a G (0, 1) and |e| < 1, one can easily check that arctan[|e| tan(7ro;/2)] G (0, na/2), and 
so na/2 < < na, < 0_ < na/2, and 0+ > 0_. 



4- a + = a_ = a G (1, 2), u + / u_ ; ^ = 

In this last case, the limiting probability density is given by the same formula (4.19). To 
find the parameters + and <f>_, we first write equations 

(U+ + U) COs(7T«/2) 
COS<y2 = :, 

yu\ + u 2 _ + 2 cos(7ra:)w + w_ 

(4.20) 

(w + — U-) sin(7ro;/2) 

sin (p — : , 

\Ju\ + uz. + 2 cos(ira)u + U- 
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which follow from Eq. (3.21). Since their solution can be represented in the same form as 
Eq. (4.17), the parameters 0+ and 0_ can also be determined from Eq. (4.18). However, 
because arctan[|e| tan(7ro;/2)] G {na/2 — n, 0), in contrast to the previous case we have 
n(a — 1) < <p + < na/2, na/2 < < n, and + < 0_. 



B. Representation of V(y) in terms of the Laplace transform 

To derive the limiting probability density in terms of the Laplace transform, in Eq. (4.8) 
we first introduce a new variable of integration r\ = (1 — r)/a and use the integral represen- 
tation r(r) = f^dze-'z*- 1 (Rer > 0) of the T function [34]. This yields 



2vr? J c _ ioo sm(nr]) 

2m J J c _ ioo \y\ sm(7r?7) \\y\ a J 



(4.21) 



with < c < min{l, 1/a}. Then, taking into account the relation 



1 f c+io ° , sin(^) „ sin ®z 

— \ dr] — r-^z = n o ( 4 - 22 ) 

2* Jc-ioo sin(7r?7) 1 + 2 cos v z + z l 

(—71 < $ < 7r) and changing in Eq. (4.21) the integration variable from z to x = z/\y\, we 
obtain the desired representation of the limiting probability density 



V(y) = - r dx q^T 5"- (4-23) 

yy! 7T J l + 2cos[(j)( y y)]x a + x 2a v ; 

Note that, although the relation (4.22) is not valid for t? = 7r, Eq. (4.23) at <p(y) = ix gives 
a correct result if 7 ? (y)| < /,( 2/ )= 7r is interpreted as the limit lim^ ^ > (y)U(2/)=7r-c- 

The limiting probability density V(y) in the form of Eq. (4.23) is useful to establish its 
general properties. In particular, directly from this representation it follows that V(y) > 
[i.e., V(y) is in fact the probability density], dV(y)/d\y\ <0(y^ 0), and ma.xV(y) = V(0). 
Moreover, because of the exponential factor in the integrand, the representation (4.23) is 
the most suitable for the numerical evaluation of V(y) at large \y\. 

For convenience of use, we write below the representations of V(y) in terms of the Laplace 
transform for all four cases considered in Sec. IV A. 
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1. a £ (0, 1), a + 7^ a_ 



According to Eqs. (4.11) and (4.23), in this case 



n Jo 1 



1 + 2 cos((() + )x a + x 



sin(0 + )a; a 



(4.24) 



with (f) + = na. Since at y — the above integral diverges, one gets V(y)\ ay ^ + o = oo. Then, 
using the standard integral [31] 



(0<|^|<7r,0<z/<2) and taking into account that at v — 1 the right-hand side of 
Eq. (4.25) equals i?/ sint?, one can easily make sure that the normalization condition for 
V(y) holds: 



The main feature of the limiting probability density in the reference case is that it is 
one-sided with V(y) = at y > if a + > a_ or at y < if a_ > ct + . These conditions show 
that V(y) is concentrated on the semi-axis where the tail index is the smallest, i.e., where 
the probability of long-distance jumps of a particle is the largest. For clarity, we note that 
the total probability of jumps in this direction, W^gn(o-) = / °° d^w[sga(a)^], can be even less 
than the total probability W^_ sgn ( cr ) = 1 — W / sg n(o-) of jumps in the opposite direction. 

The behavior of the limiting probability density at a G (0, 1) and a + ^ a_ is illustrated 
in Fig. 1. 

2. a £ (1, 2), a + ^ a-, h = 

For these conditions, the limiting probability density in terms of the Laplace transform 
reads 




(4.25) 





(4.26) 




+ H{-ay)-e-\ y \ 



a 



(4.27) 
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FIG. 1: (Color online) The limiting probability density at a+ = a = 1/2 and a_ > a+. The 
solid blue line represents the theoretical result (4.24) and the simulation results (see Sec. VI) are 
indicated by red triangles. 



where 
and 



7t(q! — 1). It can be shown with the help of Eq. (4.25) that V(y) is normalized 



V(y)\ 



cry-t+O ' 



dx 



sin 



}x 



1 + 2cos(0 + )x° + x 2a 



-1 r dz , 

na Jo 1 + 2 cos 



)z + z 2 



sin((j) + /a) 



Since 



a sin(7r/a) 

7r(a — 1) and, in accordance with Eq. (4.27), V{y)\ a y^-o = 1/a, we obtain 



(4.28) 



V(0)=V(y)\^ ± o = -. 



(4.29) 



Thus, in contrast to the previous case, the limiting probability density is two-sided and 
is bounded at the origin. One branch of V(y), left if a_ > a + or right if a + > a_, is 
purely exponential, and the other is heavy-tailed (see also Sec. V). As before, the latter 
is concentrated on the semi-axis where the tail index is the smallest. Interestingly, the 
probability J °° dyV[— sgti(a)y] = 1/a that aY (oo) < 0, i.e., the total probability defined by 
the exponential branch, is larger than 1/2. 

Figure 2 illustrates the behavior of V(y) in this case. 
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FIG. 2: (Color online) The limiting probability density at a + = a = 5/4 and a_ > a+. The solid 
line with an exponential left branch is obtained from Eq. (4.27) and triangles show the simulation 
results. 

3. a + = a_ = a £ (0, 1), u + ^ u_ 



where the parameters <p + and 0_ are given by Eq. (4.18). Since a < 1, one has V(0) = oo 
and, using again Eq. (4.25), it can be verified that V(y) is normalized. 

The comparison with the first case shows that, while the difference in the tail indexes 
a + and a_ leads to a strongly asymmetric one-sided V(y), the difference in the parameters 
u + and U- (under condition that a + = a J) results in a less asymmetric two-sided V(y). 
According to Eq. (4.30), both branches of the limiting probability density have heavy tails 
characterized by the same tail index a (see also Sec. V). 

In this case, the behavior of V{y) is illustrated in Fig. 3. 

4- a + = a_ = a G (1, 2), u + ^ n_ ; l\ = 

As in the previous section, the limiting probability density and the parameters 0+ and 
0_ are determined by Eqs. (4.30) and (4.18), respectively. The striking difference between 



From Eqs. (4.19) and (4.23) it follows that 





(4.30) 
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FIG. 3: (Color online) The limiting probability density at a+ = a_ = a = 1/2 and e = 1/3. 
The solid line and triangles represent the limiting probability density (4.30) and simulation data, 
respectively. 

the behavior of V(y) in these cases is that now V(0) < oo. To find V(0), we use Eq. (4.30), 
which together with the standard integral (4.25) yields 

sm((f>±/a) 



V(y)\ ( 



a sin(7r/aj 

Using Eq. (4.18), one obtains sin(0 + /a) = sin(0_/a), and so V(y)\ ey ^±o = V(0), where 



(4.31) 



V(0) 



1 



cos^ — arctan 



e tan — 
11 V 2 



(4.32) 



asin(7r/a) \a 

Comparing with the second case, we again observe that the difference in a + and a_ causes 
more change in the behavior of branches of the limiting probability density than the differ- 
ence in u + and w_ at a + = a_. 

The behavior of V(y) in this last case is shown in Fig. 4. 



C. Representation of V(y) in terms of the Fox H function 

Since the Fox H function is one of the most general special functions and many of its 
properties are well studied, it is also reasonable to express the limiting probability density in 
terms of this function. For this purpose, we first use the reflection formula [34] T(z)T(l—z) = 
7r/ sin(7rz) to obtain 

sin[0(y)^] _ r(^)r(i-^) 

sin(Tr^) r[0(y)^]r[l-0(y)^]- 
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FIG. 4: (Color online) The limiting probability density at a + = a_ = a = 5/4 and e = (y/2 - 
1)/(V2+1). As in the previous case, the solid line is obtained from Eq. (4.30) and triangles indicate 
the simulation results. 



Then, substituting this relation into Eq. (4.8), one gets 

1 



2iria 



c+ico r(r)r(i=^)r(i 

dr- ' " 



l-r 



) 



(4.34) 



r[^)^]r[i-0(y)^] m ' 

On the other hand, the H function can be defined in the form of a Mellin-Barnes integral 
as follows [36]: 



Tjin,n 



[a p , Ap) 
K B q ) 



Tjin,n 
H p,q 



2ni 



(d, At) 
(bx,B 1 ) 

drQ r y~ r , 



• ; (%>; Ap) 

, K B q ) 



where 



6, 



117=1^ + Bjr) ULiTil-aj- A 3 r) 



(4.35) 



(4.36) 



i / — m+ i r(i - b 3 - B f ) YiU+i r (% + V) ' 

m, n, p, q are whole numbers, < m < q, < n < p, a 3 and bj are real or complex numbers, 
Aj, Bj > 0, L is a suitable contour in the complex r-plane which separates the poles of the 
gamma functions T(bj + B 3 r) from the poles of the gamma functions T(l — a 3 - — A 3 r\ and 
the empty product is assumed to be equal to 1. Therefore, comparing Eq. (4.36) with the 
integrand in Eq. (4.34), we obtain the following representation of the limiting probability 
density through the H function: 



-H, 



a 



2,1 

2,3 



'1 — I I) (1 



4>(y) <My) ~ 

7TQ; ' Tva > 



(0,l),(l-i i),(l 



7ra ' 7tq ' 



(4.37) 
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Using this general formula, it is not difficult to find the corresponding representations for 
all four cases considered above. But here we focus on the first one for which Eq. (4.37) can 
be further simplified. Indeed, since in this case (f>(y) = (f) sga (ay) with <p + = na and 0_ = 0, 
from Eq. (4.37) and the reduction formula 



irm,n 



(ai, A), . . . , (ap_i, A p _i), (6i, B) 
(b 1 ,B 1 ),...,(b q ,B q ) 



Tj-m—l,n 
- U p-l,q-l 



(ai, A), (ap_i, A p _) 
(b 2 ,B 2 ),...,(b q ,B q ) 



(4.38) 



[m > l,p > n) we obtain 



# M „i,i 

■-"1,2 



(1-1 1), (0,1) 
(0,1), (1-1,1), (0,1). 

(1-1 1), (0,1). 



(4.39) 



Then, using the relation 



TTm,n 

p,q 



(b q ,B q ) 



X TTm,n 



y 



Tir, 
Ax P'l 



(a p + \xA p , xA p ) 
(b q + \ X B q , X B q ) 



(4.40) 



(x > 0, — oo < A < oo) with x = a an d X — l/a — 1, Eq. (4.39) can be reduced to 



(0,1) 

(0,1), (!-«,«) 



(4.41) 



#1,2 



(4.42) 



Finally, since this if function is closely related to the generalized Mittag-Lefner function 

E a Az) [36], 

(0,1) 

(0,1), (1-/3, a) 

(a, (3 > 0), for the limiting probability density in the considered case, when a G (0, 1) and 

a + ^ a_, one gets 

nv) = T0zE Oha {-\y\"). (4.43) 

\y\ 

The usefulness of this result arises from that the Mittag-Leffler function is well studied 
(see, e.g., Ref. [37] and references therein). In particular, using the series definition of this 
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function, E a ^(z) = Y^=o z n /T(an + we obtain 



■l) n \y\ 



n=0 



F[a(n + 1)]' 



(4.44) 



Remarkably, the limiting probability density (4.43) at a — 1/2 can be expressed through 
a simple complementary error function erfc(z) = (2/y/ir) f^°dxe~ x2 . To show this, we first 
note that 



E a A-\y\ a ) = -\v 



l-Q 



^i(-br)- 



(4.45) 



(This relation follows directly from the series representation of the Mittag-Leffler function.) 
Then, using the known result [37] E 1/2 ,i(-z) = e z erfc(z), Eq. (4.43) at a = 1/2 can be 
written in the form 

(4.46) 



V{y) = H(ay) 
The plot of this density function is shown in Fig. 1. 



, — e' y 'erfc(A/j?/j) 

v n \y\ 



D. Series representation of V(y) 

We complete our study of alternative forms of the limiting probability density V(y) 
by determining its series representation. This representation can be useful, for example, 
for the numerical evaluation of V(y), especially in the vicinity of small \y\. Our starting 
point is the limiting probability density written in the form of inverse Mellin transform 
V( V ) = (2^i)~ l J^dr,SM with 

SM = r(i- a ,)2!M|,,r-> (4.47) 

Sin(7T77) 

[see the first line of Eq. (4.21)]. To calculate the above integral, we close the integration 
path by a semicircle Cr of a large radius R, which lies in the right half-plane of the complex 
variable 77. If this semicircle does not cross any singularity of S(i]), then, using the Stirling 
approximation for the T function [34], it can be shown that the contribution of Cr into the 
integral over the closed contour L vanishes as R — > 00. Therefore, from the residue theorem 
(see, e.g., Ref. [38]) we obtain V(y) = — J2j Res(5', rjj), where Res(S', rjj) denotes the residue 
of S(rj) at r] = rjj, the sum is taken over all isolated singularities (in our case poles) of S(rj) 
inside the contour L, and the sign " — " accounts for the direction of L. 
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According to Eq. (4.47), the poles of S(r)) result from the first-order poles i] n = n/a 
(n > 1) of T(l — arf) and from the first-order poles r\ m = m (m > 1) of 1/ sin(7r?7). If a is 
irrational, then these sets of poles, {r] n } and {i] m }, are not intersected, and hence all poles 
of S(r]) are also first-order. However, if a is rational, then some (or all if a — 1) poles from 
the set {r) n } coincide with some (or all) poles from the set {i] m }, resulting in the appearance 
of the second-order poles of S(rj). Since the cases with irrational and rational values of a 
seem quite different, we consider them separately. 



1. Irrational values of a 

In this case, the limiting probability density is written as V(y) = — Yl^=i Res(<S', n/a) — 
Ylm=i R> es (S, m). Therefore, taking into account that T(l — ai])\ v=n / a+ ^ ~ (— l) n /[ar(n)£] 
and 1/ sin(7r?7)| r?=m+ ^ ~ (— l) m /(7r£) as £ — > 0, and using the reflection formula T(l — am) = 
7r/[r(am) sin(7ram)], we readily find 



W> a ^ r(n) sin(7rn/a) 

n=l v ' 



\y\ 



+ £. (-i)-'^H 

z — ' I am sin (ream) 

rn=l \ J \ ! 

2. Rational values of a 

Let us assume now that the tail parameter a is given by the irreducible fraction a = l/p, 
where /(> 1) and p(> 1) are natural numbers satisfying the condition / < 2p. In this case, 
the first-order poles of T(l — prj/l) with numbers n = Ik (k — 1, 2, . . .) and the first-order 
poles of 1/ sm^irr]) with numbers m = pk are merged, and thus the poles of S(rj) at r\ = pk 
become second-order. It is therefore convenient to represent the limiting probability density 
in the form 



p( y ) = -J2Res(S,pn/l) - Res(S,m) 

n=l m=l 
(n=il,2l,...) (m^p,2p,...) 

oo 

-^Res(S,pA;), (4.49) 

k=l 

where the last sum is over all second-order poles of S(k). Using the above results for the 
residues of S(rj) at the first-order poles and the asymptotic formula [34] r(l — lrj/p)\v=pkH ~ 
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(-l^p/Otr 1 - (l/p)il>(lk)]/r(lk) (f ->• O) with V(z) = rflnr(x)/rfa; being the ^ (or 

digamma) function, from Eq. (4.49) we obtain 

py^ (-l)"" 1 sm[(j)(y)pn /l] LJn _ x 
l ~[ T(n) sin(7rpn//) 



(n^,2/,...) 

+ (-l) m ~ 1 sin[0(y)m] ., ro/H 
^-^Y{lm/p) sin(7r/m/p) 



(m,y^p,2p,...) 



-^(^oos^pfc])^!'*- 1 . (4.50) 

It should be noted that if a e (0, 1) and a + 7^ a_ then the function <f)(y) is given by 
Eq. (4.10). In this case Eqs. (4.48) and (4.50) are reduced to Eq. (4.44) leading to the 
Mittag-Lefner function representation (4.43). If a G (1,2) and a + 7^ a_ then V(y) can 
also be expressed in terms of the Mittag-Lemer function. Indeed, using Eq. (4.12) and the 
relations £~ (±l) n M7™! = e ±|j/| and £~ =1 \y\ an /T(an) = E afi (\y\ a ) = \y\ a E a , a (\y\ a ), 
Eqs. (4.48) and (4.50) can easily be reduced to 



v\ 



V(y) = H(<jy)^—-\y\ a -'E a , a (\y\ a ] 

e -\y\ 

+ H(-ay) . (4.51) 

a 

V. SHORT- AND LONG-DISTANCE BEHAVIOR OF V{y) 

The short-distance behavior of the limiting probability density V(y) is completely de- 
scribed by the series representations (4.48) and (4.50). To find the long-distance behavior 
of V(y), it is convenient to use its Laplace transform representation (4.23). According to 
Watson's lemma [38], the asymptotic series expansion of V(y) at \y\ — > 00 is determined 
from the small- a; series expansion of the integrand function multiplied by e' y K Therefore, 
using the series expansion [39] 

= y(_l)"-i s in[0(y)n]^" (5.1) 

l + 2cos[(p( y y)]x a + x 2a ^ ! [Vyy! J V ; 

(\x\ < 1) and the standard integral [31] 

rdxe-l^= r ( 1 + aw) , (5.2) 

Jo \y\ 1+an 
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from Eq. (4.23) one obtains 

nv) - l - B-ir 1 Mm^^- (5.3) 

n=l |y| 

as |y| — > oo. For all cases of interest, we list below the main terms of V(y) at \y\ — > and 
\y\ — >■ oo. 

a G (0, 1), a+ 7^ a_ 

In this case, the two-valued function (f>(y) is given by Eq. (4.10), V(y)\ ay <o = 0, and thus 
Eqs. (4.44) and (5.3) lead to 

P(y)U><>~ J x, M_ a ( 5 - 4 ) 

1 (a) li/l 1 a 

(|y| — )• 0) and to 

1 1 
^(y)U>o ~ -sin(7ra)r(l + a)^— ^ (5.5) 

(|y| — > oo), respectively. 

2. a G (1, 2), a + / a_ ; ^ = 

For these conditions, the function <p(y) is determined from Eq. (4.12), V(y)\ ay <o = e~^/a, 
and so Eqs. (4.48), (4.50) and (5.3) yield 

WL^-^-^br 1 (5.6) 
ct 1 (a) 

as 1 2/ 1 — > and 

^(y)U>o ~ --sin(7ra)r(l + a)^— (5.7) 

7T 2/ 

as 1 2/ 1 — > oo. 

5. a + = a_ = a G (0, 1), -u + / u_ 

Using Eq. (4.18), it can be straightforwardly shown that 

P(V) / + Sgn(e?/)|e| , ,1 (5.8) 

2T{a)^e 2 + (1 - e 2 ) cos 2 (7ra/2) M^ Q 
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as \y\ — > and 

_ [1 + sgn(ey) |e|] sin(7ra)r(l + a) 1 , g g , 

2ny/e 2 + (1 -e 2 )cos 2 (7ra/2) |y| 1+Q 
as |y| — > oo. In contrast to the first case, the limiting probability density has both left and 
right branches characterized by the same tail index a. 



4- a+ = a- = a £ (1, 2), u + ^ U-, l\ = 
Finally, in this case we have 

nv)~V( 0) - ] + ^f} |, r - (5.10) 

2r(a)A/e 2 + (1 — e 2 ) cos 2 (7ra/2) 

x [1 + sgn(ey) |e|] sin(7ra)r(l + a) 1 ^ 

27r v /e 2 + (l-e 2 )cos 2 (7ra/2) 
as |y| ->■ oo, where P(0) is given by Eq. (4.32). Note that at e = Eqs. (5.8)-(5.11) are 
reduced to those obtained in Ref. [29] for symmetric walks. 



as \y\ — > and 



VI. NUMERICAL SIMULATION OF V{y) 

The determination of the limiting probability density V(y) by the numerical simulation 
is not a trivial problem. To understand why this is so, we first recall that V(y) is the 
probability density of the random variable Y(t) = a(t)X(t) in the limit t — > oo. In the 
simulation, however, we need to consider the behavior of the variable Y(T) = a(T)X(T) 
and the corresponding probability density 

for some finite time t = T. To be sure that this probability density approaches V(y), the 
operating time T must be large enough and, in principle, it should exceed the characteristic 
scale of waiting times. But in our case all fractional moments of the waiting-time density 
p(r) do not exist and thus there is no finite time scale of p(r). This means that for any finite 
T there is always a non- negligible survival probability V(T) = J™ drp(r) that the waiting 
time is larger than T. Therefore, the minimal value of T is restricted only by the condition 
a(T) <^ 1 which is equivalent to V(T) <^ 1. Since V(T) is a slowly varying function, it 
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decreases with increasing of T very slowly, and thus the operating time is expected to be 
very large. On the other hand, the larger is T the larger is the average number N(T) of jumps 
occurring in the time interval (0, T) [this is so because, according to [40], N(T) ~ V -1 (T')], 
and hence the larger is the computational time. Thus, the chosen value of the operating 
time T must satisfy the condition V(T) <C 1 and provide a reasonable computational time. 
In our numerical simulations, we use the following waiting-time probability density: 



vWg 
( g + r)\n 1+v (g + r) 



P(r) = T , Jf, , - (6-2) 



with v > and g > 1. The main advantage of this density is that its distribution function 
-^p( r ) = Jo dr'p{r') is calculated explicitly 

/ \ In 11 q 

F < M = 1 - wyh-y (6 - 3) 

and thus the inverse function of U — F p (r) is given by r = g^ u ^> 1/v — g. The last result 
permits us to use the inversion method [41] in accordance with which the random variables 
defined as 

r n = g^- 1/v - g, (6.4) 

where n — 1,2,... and U n are random numbers uniformly distributed in the interval [0, 1], 
have the same probability density (6.2). This provides a simple way to generate the waiting 
times. For the simulations we choose g = 2, v = 2 and T = 10 15 yielding V(T) rs 4 • 10~ 4 . 

Our choice of the jump density w(£) is limited by two conditions. First, to verify the 
theoretical results, w(£) must reproduce all possible cases considered earlier and, second, 
to simplify the generation of the jump lengths £ n , the corresponding distribution function 
F w (£) = f_ d£'w(£') must be invertible. These conditions are satisfied, for example, by the 
jump density 

fa_c_6^-/(6_-0 1+a -, e<0 
w(0 = \ (6.5) 

[a + c + b a + + /(b + + 1+a+ , e>0. 



Here, b± G (0, oo) and c + + c_ = 1 with c + and c_ being the probabilities that £ > and 
£ < 0, respectively. It can be easily shown from Eq. (6.5) that 

c_b°L-/(b- -£)<*-, £<0 
^(0 = < (6.6) 
l- c+ b a + + /(b + + O a+ , £>o, 
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and thus the jump lengths can be determined as 

[ -6_(c_/C/ n ) 1/a - +b-, U n <c_ 

^ = { (6.7) 

[& + [ c+ /(i-iy] 1/a+ -&+, c/„> c _. 

Now we are ready to describe the procedure for calculating Vt{v)- According to the 
definition, the particle starts to walk at time t — from the position X(0) = 0. At 
the first step, the waiting time T\ and the jump length £x are generated using Eqs. (6.4) 
and (6.7), respectively. If T\ < T, then the particle position becomes X(t±) = £i, and 
we can go to second step that consists in generating new random numbers r 2 and £ 2 - If 
at the nth step the condition XT=i r * — ^ * s violated, then the walk is stopped at the 
previous step, i.e., the scaled position Y(t) of the first particle at t = T is assumed to be 
Y(T) = a(T) YllZi it IX(^) = at n = 1], where the scaling function a(T) can be calculated 
using an appropriate theoretical formula. Determining Y(T) for N particles, we can evaluate 
the limiting probability density as follows: Vt(v) = Na v /N, where N& y is the number of 
particles with Y(T) e [y,y + Ay). In all our simulations we set = 10 5 and Ay = 1CT 1 ; 
the other parameters are listed below. 



1. a £ (0, 1), a + 7^ ct_ 

In this case, the limiting probability density V(y) depends only on the minimal tail index 
a [see, e.g., Eq. (4.11)]. But to determine the approximate probability density Vt(v) by 
the proposed procedure, all the parameters in Eqs. (6.2) and (6.5) must be specified. In 
addition to the parameters mentioned above, we choose a = a + = 1/2, c + = 2/3, b + — 1 
and a_ = 3/4, c_ = 1/3, = 1. With these parameters, Eq. (3.18) yields a{T) w 1.2- 10" 7 
and the simulated values of Vriy), marked by red triangles, are shown in Fig. 1. 



2. a £ (1, 2), a + / a_, h = 

Since the first moment of the probability density (6.5) is given by 

/ c+h + c ~ b ~ (a z\ 

a + — 1 a- — 1 

(a± > 1) the condition l± — implies that the parameters of w(£) satisfy the condition 
c + b + /(a + — 1) = c_6_/(a_ — 1). In particular, if a = a + = 5/4, c + = 5/22, b + — 1 
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and a_ = 37/20, c_ = 17/22, 6_ = 1, then c+b + /(a+ - 1) = c_6_/(ai_ - 1) = 10/11, 
a(T) ps 1.8 • 10~ 3 , and the results of simulation of Vt(v) are shown in Fig. 2. 

5. a + = a_ = a G (0, 1), u + ^ m_ 

According to Eq. (6.5) and the asymptotic formula (2.13), the parameters u± are ex- 
pressed through the parameters of w(£) as follows: u± = a±c±b± t . Keeping in mind that 
a + = a_ = a and u + ^ U-, i.e., c+6+ + ^ cJb a S , we choose a = 1/2, c + = 2/3, b + — 1 and 
c_ = 1/3, = 1. For these parameters, w+ = 1/3, u_ = 1/6 (i.e., e = 1/3), a(T) w 9.3 -10" 8 
and the simulated results are shown in Fig. 3. 

4- a+ = a_ = a G (1, 2), u + 7^ u_ ; /1 = 

In this case, the conditions u + 7^ w_ and li — are reduced to c + 6" 7^ c_6° and 
c + 6 + = c_5_, respectively. Choosing a = 5/4, c + = 1/5, b + = 5 and c_ = 4/5, 5_ = 5/4, 
one gets u + = 5 5 / 4 /4 « 1.87, m_ = (5/4) 5 / 4 « 1.32 [i.e., e = (v 7 ^ - 1)/(V2 + 1) ~ 0.17], 
a(T) pa 5.2 • 10~ 4 and the simulated values of Vt{v) are shown in Fig. 4. 

As is seen from these figures, the numerical results are in very good agreement with 
our theoretical predictions. It is also worth to note that the proposed numerical method 
reproduces all other theoretical results [e.g., the limiting probability density at l± 7^ and 
w(—£) = w(£)] and can easily be extended to study the CTRWs, including coupled ones, in 
higher dimensions. 

VII. CONCLUSIONS 

We have studied in detail the long-time behavior of the decoupled CTRWs characterized 
by superheavy-tailed distributions of waiting times and asymmetric heavy-tailed distribu- 
tions of jump lengths. The main attention is devoted to introducing the scaled particle 
position and deriving its limiting probability density V(y). Using the Montroll- Weiss equa- 
tion in the Fourier-Laplace space and the asymptotic properties of the waiting-time and 
jump-length distributions, we have found both the scaling function, which determines the 
scaled position, and the representation of V(y) in terms of the inverse Fourier transform. 
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It has been shown that while the scaling function depends on the parameters describing 
the asymptotic behavior of both waiting-time and jump-length distributions, the limiting 
probability density is completely characterized by the parameters of the latter distribution. 
Among these parameters, the main role plays the smallest tail index a. 

To get more information about the limiting probability density, we have derived a number 
of alternative representations of V(y). The representation of V(y) in terms of the inverse 
Mellin transform is important from a theoretical point of view (all other representations 
considered in the paper follow from this one) and permits to determine the intervals of a 
where V(y) exhibits qualitatively different behavior. We have also obtained the limiting 
probability density in terms of the Fox H function. The importance of this representation 
is that the H function is well studied and many of the special functions can be considered 
as its particular cases. In particular, we have shown that if the tail indexes of the jump 
density are different and a G (0, 1) then V(y) is expressed through the generalized Mittag- 
Leffler function. Then, using the limiting density in terms of the Laplace transform, we 
have analytically demonstrated that V(y) is non- negative, has a maximum value at the 
origin, and monotonically decreases (or equals zero) as \y\ increases. We have also derived 
the series representation of V(y) which, in the case when the tail indexes of the jump 
density are different and a G (1,2), has been used to obtain the limiting density in terms 
of the Mittag-Leffler function. Moreover, the series representation of V(y) together with its 
Laplace transform representation have permitted us to completely describe the short- and 
long-distance behavior of the limiting probability density. It has been shown, in particular, 
that the tail index of V(y) is equal to the lowest tail index of the jump probability density. 

Finally, we have developed a numerical method for calculating the limiting probability 
density. This method, which deals with the statistics of the scaled particle position at large 
times, has been applied to calculate V(y) in all cases of interest. It has been shown that the 
simulation results for the limiting probability density are in excellent agreement with our 
theoretical predictions. 
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Appendix: Fractional equation for V(y) 

Let us define the Riesz-Feller space-fractional derivative y Dl of order 7 and skewness 9 
as (see, e.g., Ref. [42]) 

HyDV'iy)} = -e^^ e ' 2 \K\y K , (A.l) 

where 7 e (0,2], |0| < min{7,2- 7 }, and F{f(y)} = f^dye^ f(y) = f K . Then, using this 
definition and Eq. (3.19), one gets 

HyD a ^V{y)} = -e-^^\K\ a V K 

= -*(k)V k . (A.2) 

Finally, taking into account that J r {S(y)} = 1 and, as it follows from Eq. (3.5), J z {V(y)} = 
[1 + the fractional equation for the limiting probability density V(y) can be written 

in the form 

y D« 2 ^V(y) = V(y)-5(y). (A.3) 
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